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DEVELOPMENT OF A NEW TECHNIQUE FOR IMAGE 
RECONSTRUCTION, ENHANCEMENT, AND VISUALIZATION 


Principal Investigator: Dr. Steven A. Orszag 

1. INTRODUCTION 

The problem of describing and reconstructing a surface, both analytically and 
visually, has been extensively addressed in this work. Various algorithms have been 
developed, tested, and compared and several conclusions have been reached. They are: 
(1) If the surface information is given on a uniform grid, the surface is best described by 
a set of contour plotting routines that exploit the GL graphics libraries both in a solid as 
well as a wire-mesh framework; (2) If the surface information is given on a non-uniform 
grid, the surface can be described by Delaunay triangulation, a sophisticated graphics 
computer, or by interpolation back to a uniform grid. Each of these methods possesses 
certain advantages over the others. In the event of insufficient information on the 
surface, two basic methods have been employed: (1) an algorithm to fill in "missing" 
information that employs Delaunay triangulation; and (2) interpolation amongst various 
grid points of known data assuming that the surface is relatively smooth. In the case of 
distinct missing information on outliers (regions of the surface possessing extremal 


properties missing from the original surface), no progress has been made. 


2. TECHNICAL BACKGROUND 


In general, a surface is described by a function, such as f(x,y,z)=const 


Alternatively, a surface can be described by z = g(x,y), with the understanding that 


g(x,y) can take on multiple values of the surface folds back on itself. A typical area of 
interest to the Air Force is a surface reconstruction of photographs taken from high 
altitude. In some portions of the photographs where high resolution is needed, certain 


blurriness associated with absence of resolution may need to be improved. Furthermore, 
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the need to remove some white noise from the picture may result in a distinct sharpening 
of the picture. Statistically, as long as the surfaces involved are very smooth, standard 
regression analysis could very successfully yield the desired results. However, the 
information that one wants to obtain are more often than not associated with some 
characteristics that are absent from the general smooth background of the images. 

There are two distinct problems in surface reconstruction: 

(1) Decomposition of the surface into elements. 

(ii) Smoothly connecting these elements by interpolation. 

The best solution to these problems is different for uniform grids than it is for non- 
uniform grids. 

For non-uniform grids, the first problem can be solved by the Voronoi-Delaunay 
construction. Each vertex is surrounded by a Voronoi cell, defined as that cell that 
contains points which are closer to this vertex than to any other. The boundaries 
between neighboring Voronoi cells are lines, equidistant from a pair of vertices. 
Connecting all such pairs by links generates a Delaunay triangulation. We exploit this 
geometrical construction by observing that each triangle (abc) involves vertices which are 
Voronoi neighbors. There is a point Y where all three Voronoi cells (a,b,c) touch. By 
constriction, all three vertices a,b,c are equidistant from Y so that these points are on the 
circle centered on Y. No other vertex can exist inside this circle by definition of the 
Voronoi cells. This property has been used in the construction of an incremental 
algorithm for Delaunay triangulation. We have developed an algorithm that accomplishes 
this triangulation recursively. Strictly speaking, the task is to add a point to an existing 
triangulation of N points. An added point exists and is contained inside one of the 
existing triangles. If we connect the new point to the vertices of the triangle in which it 
lies, we get the zeroth approximation of the triangulation of N+1 points. This 
triangulation is good away from the new point, but must be rebuilt in the local 


neighborhood of the new point. The maximum area of reconstruction is given by the 


| 


circle centered on Y, and in general, it requires only a finite number of bond flips 
(triangle reconstructions). A recursive flipping algorithm can be derived and is very 
fast. The most costly operation in this incremental Delaunay triangulation is the search 
for the triangle surrounding the new point. Straightforward scanning takes O(N) time, 


but preprocessing the data accelerates the search by the following strategy. 


One can order the vertices along any specific curve so that neighbors on the curve 
will be always close to each other on the plane. (This is a version of the traveling 
salesmen problem.) Such an ordering of initial points takes O(N log N) time. This 
ordering is needed to be performed only once, after which various progressive surfaces 
can be drawn using the same coordinate system. Degenerate cases, when the points are 
too close to each other or the (N+1)st point is too close to the side of the Nth 


triangulation, are treated separately to prevent numerical instabilities. 


The second approach is based on the graphics implemented in the graphics 
terminals of Evans and Sutherland. The points are loaded on the system and are 
displayed simultaneously with x,y, z coordinates all belonging to the same surface. The 
surface can be rotated around each of the three axes and cutting planes have been 
introduced to isolate specific projections of the surface on a plane of choice. On a plane 
of choice, curves have been interpolated using optimal bicubic splines. The points are 
added on the splines and the files of the individualized planes get re-assembled. A 
perpendicular plane to the original set of parallel planes is being introduced and the added 
points in each of the previous planes is compared to the initial set of points on this cut 
plane. A least-squares estimator is performed to decide on admissibility or 
inadmissibility of the added points. If the set of points added to the global grid is 
admissible, a point-by-point elimination process for the added points is invoked until the 


set of added points becomes admissible. The set of points that are then admissible per 


perpendicular plane are then added to the global grid and the process may be repeated 
until no more points are admissible. At that stage, the entire file is re-loaded on the 
graphics system and a surface reappears. The surface can be described by the actual 
collection of dots, wire-mesh, or by standard graphics shading routines on the graphics 
processors. We are now in the process of adapting these algorithms to the SGI GL 


libraries. 


The third case for the non-uniform grid can be handled in two different ways. 
First, a uniform grid is imposed in the x-y plane. At each grid point, a characteristic 
radius is introduced. Each point in the x-y plane, xj, yj is taken into account. The 
number of points p in the enclosed domain is recorded and a pth Lagrange interpolation 
scheme is introduced to compute the value z associate with the said grid point. When all 
the uniform grid point values have been obtained, the interpolation scheme is checked 
against each existing original point. In other words, each original point is attempted to 
be recomputed to check for accuracy. This is done with variable orders in order for 
outliers estimates to be rejected. In the second case, at any vertex in the Delaunay 
triangulation, an arbitrary number of triangles can meet and the property of smoothness 
at ihe vertex requires a large number of local parameters. We have the following 
construction that satisfies all the required smootaness conditions. The Gauss formula 
for planes of the corresponding triangles meeting at a vertex has been used. The smooth 
function so obtained passes through all its neighbors. This function is usually a sum of 
about six terms. Three of such functions are connected inside each triangle by a similar 
Gauss interpolation resulting in an explicit local set of rational functions within each 
triangle matching their neighbors with any prescribed smoothness. The addition of the 
interpolated points to the original triangulated surface is composed by the set of points at 
the center of the Delaunay circles, i.e. vertices of the Voronoi graph which are always 


located in the voids of the Delaunay triangulation. Adding these circle centers does not 


require a search and is numerically stable. It roughly doubles the number of points 
omitting the points whose Delaunay circle radius is smaller than a given criterion to avoid 
condensation of points. After a few steps, the plane is completely covered and the 


process is complete. 


The issue of describing a surface on a uniform grid is different in nature. An aerial 
photography or topographic description of a surface can in principle be described on a 
uniform grid since the data gathering can be digitized according to the data gatherer. At 
this juncture, a missing observation on the uniform grid is interpolated using standard 
Lagrange interpolations and data contained within noise can be filtered out by invoking 
standard FFT techniques. In particular, a combination of low and high-pass filters can 
isolate specific structures by sharpening the domain of resolution. It is believed that this 


technique may have the potential of image enhancement. 


One of the major goals in this study is the ability to characterize an irregular 
su.face, in particular, non-uniform surfaces obtained either from exterior sets of 
observations or from solutions of partial differential equations of moving fronts or 
equipotential surfaces. For example, the system of detecting irregular conductivity of an 
internal human tissue is dependent on surface potential and conductivity. A solution of 
the potential problem inside the bulk with isolation of equipotential surfaces is hoped to 
improve visualization of internal irregular tissues. It is for such surfaces that the non- 


uniform grid methodologies developed here will be useful. 


3. IMPLEMENTATION 
We have developed the algorithms specified in Tasks 1-3 of the original proposal. 
These are exemplified by the attached source code in Appendix I. For Task 4, we have 


developed a menu-driven user-based interface that is displayed with all its functions in 
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the attached plots. The use of this interface and code is elementary and requires minimal 
learning time. The code interface is implemented by use of a mouse within the graphics 
window. For Task 5, we provide various sets of examples of non-uniform as well as 
uniform grid systems. The uniform grid graphics visualization displays a system 
borrowed from microelectronics of contact hole printing when the surface of the 
electromagnetic energy intensity is displayed above a given surface. Multi-resolution is 
used in these calculations to address the requirements of circuit designers. Furthermore, 
we attach two different circuit designs (see Figs. 1-2) obtained for uniform grid 
implementation of two separate logic patterns for wiring and memory triggers. In both 
cases, GEOFORM gave a reasonably satisfactory description of the circuit and its 


associated resolution. 


For the non-uniform grid systems, we have tested this method on geologic data, 
Martian topography, vortex sheets in hydrodynamics and also on a sculpture 
reconstructed from a few hundred points. Sometimes smoothing was needed, i.e., the 
points were iteratively moved towaids the local second-order surfaces, least-square fitted 
to the nearest neighbors. Since this fit was local, it did not erase the small details at the 
surface of larger structures. Everything was smoothed at its own scale. This algorithm 
is Surprisingly efficient. The Mars topography reconstruction required only 2% of 
randomly chosen data points. The sculpture was reconstructed from just 600 points, and 
the roll-up of a vortex sheet was simulated using only 1000 points. In Figure 3, we 
display the Delaunay triangulation of the initial mask. The reconstructed surface is 
displayed next to this triangulation and the internal variations in the density of the 
triangulation manifests in different shading to create the necessary visualization of the 
actual original mask. The other two segments of Figure 3 involve the reconstruction of 
very sharp mountain peaks on the Asian continent where the triangulation density 


indicated initially a sharp gradient that led to the reconstruction of the Himalayan peaks. 


‘ a = - * 
ae : , sos ant ae 5° * : . < - * 
“ + . % a 5 3 . 5 ‘ x 4 . » 
y © : re wwe ee on 7 & (ya 
‘s = * Boy oe eS PANE) alae 
vy RR Se econ oe at eee z 
aes thee oo Ss ‘ 
me ates ge Fae Ng 8 say 3 
Cen ee ee ‘¢ J 
er 8% ern ~ ears . 
a s (5946 “ xe SY wy s 
BU Se eg ean See eos Be 
Se OME & ape ‘ oe ahr ; 205 
a AME cee ; Pst: , eng ekogh $ siege Se mid 
. Ps oe es ‘ e : : i © he BAR? es as * 
wig Ste tet S ar -, 
: Oe + fat 
© Ty i 
; 4 
Sais ce Ss j eee 
Sas ee é 
es Ze PORE LOR ie ERS 
is . tee ee, 
oS fe We He ee Cw "ha 
. eat . A” She WN Tae go My Se 
often ys oy Re Sere ee oe a 
Cela) Woke s Sc ON 3 
SS nah & % A BH HM ee 
ah i¥ gh EN SS fod B. £: ay 295 y wht - tet yt + 
warts “oe e es r Cee st atta x 
set, Ww RS at BRS TS rs PO of 
Seta, Bae See Sa ty s wd eta wy te Shy 
oh, Me BS pe ae ‘ oe q 
pe ee ere Sok US Sd Hees FS : 
et ge Ee oa td At yar se 
‘ a Some S aa he les 
ee: a i a Sos S eye Mow alk gs OS 
Borge ie cebu ee S eh * 
a or oye ” ay o . 
q rer ner 5 ety ye Me 
$ " “eg : 4 : : Me 
3 Or. ane at poe Se Ago : % "1% 
te ea gy ahha, me VE Suc te wea? : é Oe : : g Sr BEC RGY M4 Bee wh 
Pea has 4 ASS $3 ia ; nt 
wl Th eNDy Jes Sp ee Pent : ae g 
Pe is Aw oy ‘ > : fe te > 
> Ti eT Ate, tet 3 $8 oe | 
‘  g Belg Seg, te : 4 
“sec hA wy hates v* ¢ oa 
cg PA 1 Bae f 
ET beer ol, Bede os Ee ; Qs 
i, 3 Re Me Ae” ag" {o> 
Ne fe ous ee & ee soe) ow 
oes mo ae ts + 
ae ‘ 
Cee em Be eae b wi fs et 
Bini reg » ; i 
oe eee be S g, % A, 086 
Ae Bho oop Re 
tee i 
ey Rhee Bae 4 
pork sg aa Ss S 
ae be 4 , pat pi. | 
are te 
tA BAL ee 4 
Fup rhe, * RQ. 
ORK . 2 
ya 1 hi re, 
vipat 5 es aya tie 
SSS ot ah Be. Fate rr e 
ekg bys Aap tew! Shy 
ie! LE ees 
ta 2 awl see PAM BO SF 
is euch ote tar Cte Os 
hI Bs Bees aS 
ca aac ra meige Pian fa 
x 4 * 
A ARRTS Bees Cosel. & 
‘ f ae 
Anis ‘ Ore a 
ve a ds. 3 z ’ 2 > ¢ 3 : : : ur bs anes 
vee 7 3% < 2 : RA ‘ Me SOR + vi 
ce he Bent. oe : ee Ea BS oe i) 
.e : a Roe : BEES Roy S94 a 
oe ee tt oat 
+, 88, Ke wt 
iS eese “Ree e3h ° 
Es * 


RNA en, 
= 


fay 


Re Nae 


aS AP ged x 

cae eas Sie tty 

38 ES &, 

% EASBTS 
s%, 


x 
SN 
se 


a 


fo 


> 
x 
ser cud 
at 
nt 
ae x 
A eee 
ad 
kas %> 
ET < 
ae 
aa 
q 


3 
x 
yt 
< 
an 
Wes 


Se 


ok 
; 
wee 


ESD 
2 Ns 


B 

i. 

oy 
SoM oo 


SN 
+ ot 


* 
eee 
a8 


4 
yee 

a 

Ne 
ae 


ee 


ae 


ate ty} 
set das 


oe phi asdcs abou seer 


3 


Y 
MD. 


rs 
x 
eS) 


: 3 De barack ; 3 EB asuesinsys Rae stg 
ee Aes fee Oe Bepeaaotan i HES x s hptas 
s : fe fd 


age 
% 


Ae 


S 
6 


Ae 
“EBA 
re , 
: PRS 
a EERE 
¥ EPs SAE 2 2 eo 
ce . oan «4 ‘ * 
< * 


1s Moe ‘ ‘ 
PE ES, Sul Ene 


peut Foxe © eer toe hoger, Abe 
2 ie Be : ioe 


Gees NS eh eee See 
< Pe oe nag 
; ZS 
+4 


SPRL aae SSNS Kg) PRET 
ERAS Ss ‘ : 8 ey FS es e Sean 3 Sa Rees 
BS Se: eek aR! 
; : ey 


es 


: 
EONS 


2 
x 


ay . : 2, bs 
WE, Ewa. Re ; or a Sia Ge 
x 2 


lg Es see ia 
son Ep secre at: 


SOOSRESE GE 

Span 

Sian 620 
es 


% 


ica 


appli 


3: Examples of GEOFORM 
tions 


Higher intensities (higher heights of the mountain peaks) are illustrated in red to 
emphasize the reconstructed of the missing peak. This technique could be very useful to 
reconstruct surfaces obtained from high altitude with blurred resolution at the extremes. 
The fourth quadrant of Figure 3 illustrates a reconstructed sphere with cavities inside. It 
could be particularly important to use this technique in the simulations of porous media 
where specific flow channels can be obtained to determine porosity of new materials. 
The dynamics of pore sizes is under current investigation. 

In Figures 4-5, we display graphically the detailed mechanism of edge detection in 
the GEOFORM code. It is menu controlled and specific points associated with the edge 
of a given surface are determined locally or globally. In particular, a zoom window 
provides the ability to zoom in on an edge of interest. The particular example given is 
that of refining the edge of a ‘rabbit' shaped region. It is apparent that even rough 
resolution suffices to solve the problem. A second graphical option associated with a 
collection of points is the imposition of the interior Delaunay triangulation embedded in 


these points. 


To explain how we do this, let us consider a collection of points which lie along 
lines. The points may be given in any order. Our objective is to divide the set of points 
into subsets, such that all the points in one subset lie along a line. In other words we 
want to reduce our initial set of points to a set of clusters of close points. We then 
compress these clusters to "skeletons". These skeletons are represented as graphs with 
logical links between points. After this compression, one can either compare graphs with 
a lookup table of standard elements or simply transmit the graph. At each link of the 
graph one can store the number and mean dispersion of the compressed nearest points, 
so that one can return to the thick lines after transmission by randomly distributing points 
alongthe links. As compared to standard methods, this one seems to be especially well 


suited for graphic objects such as handwriting, fingerprints, space photos or maps. 
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boundary of 'rabbit' 
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Figure 5: detection: refined 
boundary of rabbit using GEOFORM 


We use the following method. First, we construct the Delaunay triangulation of the 
set of given points, and then we begin a transformation of the triangulation which will 
extract the points lying along lines. We have developed a set of functions described 
below that we apply to our triangulation. The function DELETE MAX EDGE deletes the 
longest edge of the Delaunay triangulation. The function DELETE MAX EDGE IN 
TRIANGLE deletes the longest edge, which does not "disconnect" the triangulation, i.e. 
we delete the edge only if each of its end points has at least three neighbors. The function 
does nothing if such an edge does not exist. The function DELETE MIN EDGE deletes 
the shortest edgeof the Delaunay triangulation. The function MERGE MIN EDGE merges 
the end points of the shortest edge into one point. This new point is placed at the center 
of mass of the two old points taking into account the number of points which were 
merged into the points before. The function MERGE MIN TRIANGLE merges the three 
vertices of the smallest triangle (by area) into one point. The new point is placed at the 
center of mass of the three vertices, taking into account the number of points which were 
merged into the vertices before. With these functions in hand, it remains to determine 
the order of applying them to the triangulation. To explore this question, we developed a 
tool implemented on IRIS workstations which allows us to choose which function to 
apply to the traingulation at any time to achieve the best line extraction. The tool collects 
the information about all steps and prints out the statistics on the working process. This 


helps to determine the best strategy for applying our functions. 


The reconstruction of Martian surfaces is displayed in Figures 6-11. These illustrate 
extremes obtained from higher Delaunay triangulation densities with missing peaks. A 
solid surface with proper shading is reconstructed to visualize the actual surface obtained 


from discrete pieces of information. 
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Figure 6: Menu-driven wire-mes 
triangulation and refinement 
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Figure 7: Menu-based control 
environment of GEOFORM 
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Figure 8 Menu-based topo-mapping 
environment of GEOFORM 
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Figure 10: Menu-based zoom-in 
environment of GEOFORM 


Figure 11: Menu-based lighting 
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Last, but not least, we describe some of the features associated with the menu- 
driven graphics system of GEOFORM. In the GEOPLANE panel, we have a button in 
which the file name to be loaded is inserted. Then, triangulation takes place upon 
clicking on the triangulation button. After the determination of missing points takes 
place, the interpolation button gets clicked. Then the data editor is triggered and the 
display’ provides the surface. The button labeled ‘option’ allows either normalization, 


smoothing by killing peaks, or a logarithmic scale contraction. 


Another panel, the control panel, allows changes of light intensities, changes of 
colors, changes of individual rotation angles, and allows saving of files under the button 
‘Save Picture’. The color editor selects the bit of the graphics involved and allows the 
insertion of shading under the control panel instructions. These panels, as well as their 
sub-panels provide complete graphics control of the surface to be displayed either as a 
wire mesh or as a shaded solid surface. 

The remarkable flexibility of displaying results is best illustrated in the 
hydrodynamic pictures provided below. In Figures 12-15, we apply GEOFORM to the 
visualization of “ame fronts in turbulent flows. In Figures 16-20, we apply GEOFORM 
to the visualization of random of multi-frequency Rayleigh-Taylor instabilities occurring 
at density interfaces. 

In conclusion, we have developed new algorithms for data preprocessing and 
organization, triangulation and interpolation. Numerous tests verify the efficiency and 
numerical stability of these new algorithms. In particular, the surfaces are glued together 
from small elements, depending on a finite number of neighboring points, but smoothly 
connected to each other. A listing of the GEOFORM family of codes is given in the 


Appendix. 


Figure 12: Propagation of flame 
front in turbulent channei flow 
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Figure 13: Propagation of flame 
front (late time) 
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Figure 15: Head-on GE FORM view 
of flame in channel 
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Figure 16: Three-dimensional view of 
random Rayleigh-Taylor instability 
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Figure 18: Three-dimensional single 
frequency Rayleigh-Taylor instability 
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Figure 19: Early 
frequency Rayleigh: 


20: Late time 3D two- 
frequency Rayleigh-Taylor instability 


LISTING OF GEOFORM AND RELATED PROGRAMS 
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